/**
 * @file   Vector.h
 * @author Heyu Wang <hywang@walker-9.local>
 * @date   Tue Oct 27 16:08:10 2020
 * 
 * @brief A simple implement of vector space.  
 * 
 * 
 */

#ifndef __CRAZYFISH_VECTOR__
#define __CRAZYFISH_VECTOR__

#include <cmath>
#include <iostream>
#include <iomanip>


/**
 * Vector structure.
 * 
 */
template <class T,int _dim>
class Vector
{
private:
    T coord[_dim];		/**< The coordinate of a vector. */
    enum {dim = _dim};			/**< The dimension. */
public:
    /** 
     * The default constructor.
     * Use the default value of template @p T to construct the vector.
     */
    Vector();
    /** 
     * Constructor.
     * Use {coord1,coord2,...} to construct the vector.
     */
    Vector(std::initializer_list< T > l);
    /** 
     * Constructor. Provide the coordinate of the vector by a double array, 
     * and the parameters @_start and @_end are the first and the last address 
     * of the coordinate in the given array.   
     * 
     * @param _start The address of the first component.
     * @param _end The address of the last component.
     */
    Vector(T* _start, T* _end);
    /** 
     * The destructor.
     * 
     */
    //~Vector(){};//delete [] coord;};
    /** 
     * Make a mirror image of the vector.
     * The mirror image is used to avoid data loss.
     * 
     * @return The mirror vector.
     */
    Vector<T,_dim> mirror() const ;
    /** 
     * Print out the coordinate of the vector.
     * 
     */
    void print() const;
    /** 
     * Give the dimension of the vector.
     * 
     * 
     * @return The dimension.
     */
    int get_dim() const ;
    /** 
     * Give the @p _idx -th coordinate component value.
     * 
     * @param _idx The index.
     * 
     * @return The coordinate value.
     */
    const T& get_coord(int _idx) const ;

    /** 
     * Set the @p _idx -th coordinate by value @p _val .
     * 
     * @param _idx 
     * @param _val 
     */
    void set_coord(int _idx, T _val) ;
    /** 
     * Give the mathematical @p _idx -th coordinate value.
     * 
     * 
     * @return The value.
     */
    const T& operator()(int _idx) const ;
    /** 
     * Set the vector by @p _vec .
     * 
     * @param _vec The object vector for copy.
     */
    void operator=(const Vector<T,_dim> & _vec) ;
    /** 
     * Transfer to the opposite vector.
     * 
     * 
     * @return The opposite vector.
     */
    Vector<T,_dim> operator-() const ;
    /** 
     * Vector addition.
     * 
     * @param _vec The vector added to.
     * 
     * @return The sum vector.
     */
    Vector<T,_dim> operator+(const Vector<T,_dim> & _vec) const ;
    /** 
     * Vector subtraction.
     * 
     * @param _vec The vector substracted to.
     * 
     * @return The difference vector.
     */
    Vector<T,_dim> operator-(const Vector<T,_dim>& _vec) const ;
    /** 
     * Evaluate the inner product of this vector and the vector @p _vec .
     * 
     * @param _vec The vector evaluated with.
     * 
     * @return The inner product.
     */
    T dot(const Vector<T,_dim>& _vec) const ;
    /** 
     * Scalar multiplication.
     * 
     * @param c The scalar.
     * 
     * @return The product vector.
     */
    template <class T1>
    Vector<T,_dim> operator*(T1 c) const ;
    /** 
     * Scalar multiplication.
     * 
     * @param c The scalar, on the left side.
     * @param _vec The vector multipled to.
     * 
     * @return The product vector.
     */
    template <class T1,class T2,int _dim1>
    friend Vector<T1,_dim1> operator*(T2 c,const Vector<T1,_dim1> _vec);
    /** 
     * Evaluate the L1 norm of the vec.
     * 
     * @param _vec 
     * 
     * @return The norm value
     */
    template <class T1,int _dim1>
    friend T1 L1norm(const Vector<T1,_dim1> _vec);
    
    /** 
     * Evaluate the L2 norm of the vector.
     * 
     * 
     * @return The norm value.
     */
    template <int _dim1>
    friend double L2norm(const Vector<double,_dim1> _vec);
    /** 
     * Print @p _val in a given way.
     * 
     * @param _val 
     */
    template <class T1>
    friend void myprint(T1 _val);


    template <class T1>
    friend T myplus(T _vec1,T _vec2);
};
/** 
 * Transpose the matrix
 * 
 * @param mx
 * 
 * @return The transposed matrix
 */
template<class T,int m,int n>
Vector<Vector<T,m>,n> trans(Vector<Vector<T,n>,m> mx);
/** 
 * Calculate mx1*mx2
 * 
 * @param mx1
 * @param mx2
 * @return mx1*mx2
 */
template<class T,int m,int s,int n>
Vector<Vector<T,n>,m> operator*(Vector<Vector<T,s>,m> mx1,Vector<Vector<T,n>,s> mx2);

template <class T,int n>
Vector<Vector<T,n-1>,n-1> del(Vector<Vector<T,n>,n> mx,int i,int j);
/** 
 * @param mx 
 * 
 * @return det of mx
 */
template <class T,int n>
T det(Vector<Vector<T,n>,n> mx);
/** 
 * Gauss elimination 
 * 
 * @param a coefficient matrix
 * @param b constant term vector
 */
template <class T,int n>
void GElimination(Vector<Vector<T,n>,n> a,Vector<T,n> b);

/** 
 * 
 * @param A Gauss-Seidel迭代法收敛且对角线元不为0  
 * @param b
 * @param epsilon  
 */
template <class T,int n>
void GSMethod(Vector<Vector<T,n>,n> A,Vector<T,n> b,T epsilon);


////friend

template <class T,class T1,int _dim>
Vector<T,_dim> operator*(T1 c,const Vector<T,_dim> _vec) {
    return _vec * c;
}

template <class T,int dim>
T L1norm(const Vector<T,dim> _vec) {
    T result = 0;
    for (int i = 0; i < dim; ++i) {
	result += (_vec.get_coord(i) > 0) ? _vec.get_coord(i) : -1 * _vec.get_coord(i);
    }
    return result;
};

template <int _dim>
double L2norm(const Vector<double,_dim> _vec){
    double result = 0;
    for (int i = 0; i < _dim; i++)
    {
	result += double(_vec.get_coord(i)*_vec.get_coord(i));
    }
    result = sqrt(result);
    return result;
};

//myprint
// TODO: print complex number
template <class T>
void myprint(T _val){
    std::cout << "Not defined yet" << "\n";



};
template <>
void myprint(double _val) {
    std::cout << std::setw(8) << std::setprecision(12)
	      << _val << "  ";
    return;
};
template <class T,int _dim>
void myprint(Vector<T,_dim> _val) {
    for (int i = 0 ; i < _dim ; ++i) {
	myprint(_val.get_coord(i));
    }
    std::cout << "\n";
};


//plus
template <class T>
T myplus(T val1,T val2){};

template <class T,int _dim>
Vector<T,_dim> myplus(Vector<T,_dim> _vec1,Vector<T,_dim> _vec2) {
    Vector<T,_dim> result;
    for (int i = 0; i < _dim; ++i) {
	result.set_coord(i,myplus(_vec1.get_coord(i),_vec2.get_coord(i)));
    }
    return result;
};

template <>
double myplus(double _val1 , double _val2) {
    double val = _val1 + _val2;
    if (val == double (int (val))) {
	return val;
    }
    double sign = 1;
    if (val < 0) {
	val = -val;
	sign = -sign;
    }
    double bval = (abs(_val1) > abs(_val2))? abs(_val1) : abs(_val2);
    double decimal = bval;
    while (decimal > 1)	--decimal;
    bval = bval - decimal;
    
    
    double _digit = 15 - log10( (bval < 1)? 1 : bval );
    int digit = static_cast<int>(_digit);
    long long tmp = static_cast<long long> ((val + 0.5/pow(10,digit))*pow(10,digit));
    val = static_cast<double> (tmp)/pow(10,digit);
    
    return sign * val;

};
////



template <class T,int _dim>
Vector<T,_dim>::Vector()
{
    for (int i = 0; i < dim; ++i) {
	coord[i] = T();
    }
};

template <class T,int _dim>
Vector<T,_dim>::Vector(std::initializer_list< T > l){
    auto j = l.begin();
    for (int i = 0; i < dim; ++i) {
	coord[i] = *j++;
    }
 };

template <class T,int _dim>
Vector<T,_dim>::Vector(T* _start, T* _end)
{
    if (_end - _start != _dim ) {
	std::cout << "mismatch" << "\n";
	return;
    }
    T* temp = _start;
    T* p = coord;
    do
    {
        *p++ = *temp++;
    } while (temp != _end);    
}

template <class T,int _dim>
Vector<T,_dim> Vector<T,_dim>::mirror() const {
    Vector<T,_dim> mirror;
    for (int i = 0; i < dim; ++i)
	mirror.set_coord(i,coord[i]);
    return mirror;
};

template <class T,int _dim>
void Vector<T,_dim>::print() const {
    std::cout << "\n";
    for (int i = 0; i < dim; ++i)
	myprint(coord[i]);
    std::cout << "\b\b\n";

};

template <class T,int _dim>
int Vector<T,_dim>::get_dim() const{
    return dim;
};

template <class T,int _dim>
const T& Vector<T,_dim>::get_coord(int _idx) const {
    return *(coord + _idx);
};

template <class T,int _dim>
void Vector<T,_dim>::set_coord(int _idx, T _val) {
    *(coord + _idx) = _val ;
};

template <class T,int _dim>
const T& Vector<T,_dim>::operator()(int _idx) const {
    if(_idx < 1){
	std::cout << "use math coordinate in operator()" << "\n";
	return *(coord);
    };
    return *(coord + _idx - 1);
};

template <class T,int _dim>
void Vector<T,_dim>::operator=(const Vector<T,_dim> & _vec) {
    for (int i = 0; i < dim; i++)
    {
	coord[i] = _vec.get_coord(i);
    }
    return;
};

template <class T,int _dim>
Vector<T,_dim> Vector<T,_dim>::operator-() const {
    Vector<T,_dim>  result;
    for (int i = 0; i < dim; i++)
    {
	result.set_coord( i , -coord[i] );
    }
    return result;
};

template <class T,int _dim>
Vector<T,_dim> Vector<T,_dim>::operator+(const Vector<T,_dim> & _vec) const {
    Vector<T,_dim>  result;
    for (int i = 0; i < dim ; ++i) {
	result.set_coord(i,myplus(coord[i],_vec.get_coord(i)));
    }
    
    return result;
};

template <class T,int _dim>
Vector<T,_dim> Vector<T,_dim>::operator-(const Vector<T,_dim>& _vec) const {
    return operator+(-_vec);
};

template <class T,int _dim>
T Vector<T,_dim>::dot(const Vector<T,_dim>& _vec) const {
    T result = 0;
    for (int i = 0; i < dim; i++)
    {
	result += coord[i]*_vec.get_coord(i);
    }
    return result;
};

template <class T,int _dim>
template<class T1>
Vector<T,_dim> Vector<T,_dim>::operator*(T1 c) const {
    Vector<T,_dim>  result;
    for (int i = 0; i < dim; i++)
	result.set_coord(i,coord[i]*c);
    return result;
};



//






/** 
 * Set structure.
 * 
 */
template <int dim>
class Set
{
private:
    Vector<double,dim>* point;		/**< The given points */
    int n;			/**< The number of points */
public:
    /** 
     * The default constructor.
     * 
     */
    Set();
    /** 
     * Constructor. Provide the points by a Vector array, 
     * and the parameters @p _start and @p _end are the first and the last address 
     * of the points in the given array. 
     * 
     * @param _start 
     * @param _end 
     */
    Set(Vector<double,dim>* _start,Vector<double,dim>* _end);
    /** 
     * The destructor.
     * 
     */
    ~Set() {};
    /** 
     * Print out the given points.
     * 
     */
    void print() const ;

    /** 
     * Use Lagrange Interpolation Polynomial to caculate the interpolation.
     * Exist Runge Phenomenon
     * 
     * @return The interpolation
     */
    double LagrangeInPoly(double x) const ;
    /** 
     * Use Simpson Integration to caculate the definite integration of the points.
     * 
     * 
     * @return The definite integral value
     */
    double Simpson() const ;
};

template <int dim>
Set<dim>::Set(Vector<double,dim>* _start,Vector<double,dim>* _end){
    Vector<double,dim>* temp = _start;
    n = _end - _start;
    point = new Vector<double,dim>[n];
    Vector<double,dim>* p = point;
    do
    {
	*p++ = *temp++;
    } while (temp != _end);
};

template <int dim>
void Set<dim>::print() const {
    for (int i = 0; i < n; ++i) {
	std::cout << "\n" << "(" << point[i].get_coord(0)
		  << " " << point[i].get_coord(1)
		  << ")"<< "\n";
    }
};

template <>
double Set<2>::LagrangeInPoly(double x) const {
    if (n < 2) {
	std::cout << "two point at least" << "\n";
	return 0;
    }
    double sum = 0;
    for (int j = 0; j < n; ++j) {
	double fj = 1;
	for (int i = 0; i < j; ++i) {
	    fj *= (x - point[i].get_coord(0))/(point[j].get_coord(0) - point[i].get_coord(0));
	}
	for (int i = j + 1; i < n; ++i) {
	    fj *= (x - point[i].get_coord(0))/(point[j].get_coord(0) - point[i].get_coord(0));
	}
	sum += point[j].get_coord(1)*fj;
    }
    return sum;
};

template <>
double Set<2>::Simpson() const {
    if (n < 2) {
	std::cout << "two point at least" << "\n";
	return 0;
    }
    double sum = point[0].get_coord(1) + point[n-1].get_coord(1);
    for (int i = 1; i < n - 1; ++i) {
	if (i % 2) {
	    sum += 4 * point[i].get_coord(1);
	}
	else {
	    sum += 2 * point[i].get_coord(1);
	}
    }
    return (point[n-1].get_coord(0) - point[0].get_coord(0)) / 3 / (n-1) * sum ;
};



/** 
 * Matrix
 * 
 */

template <class T,int m,int n>
Vector<Vector<T,m>,n> trans(Vector<Vector<T,n>,m> mx) {
    Vector<Vector<T,m>,n> result;
    for (int i = 0; i < n; ++i) {
	Vector<T,m> tmp;
	for (int j = 0; j < m; ++j) {
	    tmp.set_coord(j,mx.get_coord(j).get_coord(i));
	}
	result.set_coord(i,tmp);
    }
    return result;
};

template <class T,int m,int s,int n>
Vector<Vector<T,n>,m> operator*(Vector<Vector<T,s>,m> mx1,Vector<Vector<T,n>,s> mx2) {
    Vector<Vector<T,n>,m> result;
    for (int i = 0; i < m; ++i) {
	Vector<T,n> tmp;
	for (int j = 0; j < n; ++j) {
	    T d = 0;
	    for (int k = 0; k < s; ++k) {
		d += mx1(i+1)(k+1)*mx2(k+1)(j+1);
	    }
	    tmp.set_coord(j,d);
	}
	result.set_coord(i,tmp);
    }

    return result;
};

template <class T,int n>
Vector<Vector<T,n-1>,n-1> del(Vector<Vector<T,n>,n> mx,int i,int j) {
    Vector<Vector<T,n-1>,n-1> result;
    if (i > n || j > n || i == 0 || j == 0) {
	std::cout << "Wrong del" << "\n";
	return result;
    }
    for (int x = 0; x < i - 1; ++x) {
	Vector<T,n-1> tmp;
	for (int y = 0; y < j - 1; ++y) {
	    tmp.set_coord(y,mx(x+1)(y+1));
	}
	for (int y = j ; y < n; ++y) {
	    tmp.set_coord(y - 1,mx(x+1)(y+1));
	}
	result.set_coord(x,tmp);
    }
    
    for (int x = i ; x < n; ++x) {
	Vector<T,n-1> tmp;
	for (int y = 0; y < j - 1; ++y) {
	    tmp.set_coord(y,mx(x+1)(y+1));
	}	
	for (int y = j ; y < n; ++y) {
	    tmp.set_coord(y - 1,mx(x+1)(y+1));
	}
	result.set_coord(x - 1,tmp);
    }

    return result;
};

template <class T,int n>
T det(Vector<Vector<T,n>,n> mx) {
    T result = 0;
    for (int i = 1; i <= n; ++i) {
	if (i % 2) result += mx(1)(i)*det(del(mx,1,i));
	else result += -mx(1)(i)*det(del(mx,1,i));
    }
    return result;
};

template <class T>
T det(Vector<Vector<T,2>,2> mx) {
    return mx(1)(1)*mx(2)(2) - mx(1)(2)*mx(2)(1);
};

template <int n>
double det(Vector<Vector<double,n>,n> mx) {
    /// check    
    double result = 1;
    double tmp;
    int sign = 1;
    for (int i = 1; i < n; ++i) {
	if (mx(i)(i) == 0) {
	    int j = i + 1;
	    for (; j <= n && mx(j)(i) == 0; ++j);
	    if (j == n + 1) {
		return 0;
	    }
	    else {
		Vector<double,n> t = mx(i);
		mx.set_coord(i-1,mx(j));
		mx.set_coord(j-1,t);
		sign *= -1;
	    }
	}
	for (int j = i + 1;j <=n ; ++j) {
	    tmp = mx(j)(i)/mx(i)(i);
	    mx.set_coord(j - 1,mx(j) - tmp * mx(i));
	}
    }
    
    for (int i = 1; i <= n; ++i) result *= mx(i)(i);
    return sign * result;
};

template <class T,int n>
void GElimination(Vector<Vector<T,n>,n> a,Vector<T,n> b) {
    /// check and so on
    Vector<T,n> x;
    T tmp;
    for (int i = 1; i < n; ++i) {
	if (a(i)(i) == 0) {
	    int j = i + 1;
	    for (; j <= n && a(j)(i) == 0; ++j);
	    if (j == n + 1) {
		std::cout << "不存在唯一解" << "\n";
		return;
	    }
	    else {
		Vector<T,n> t = a(i);
		a.set_coord(i-1,a(j));
		a.set_coord(j-1,t);
		tmp = b(i);
		b.set_coord(i-1,b(j));
		b.set_coord(j-1,tmp);
	    }
	}
	for (int j = i + 1;j <=n ; ++j) {
	    tmp = a(j)(i)/a(i)(i);
	    a.set_coord(j - 1,a(j) - tmp * a(i));
	    b.set_coord(j - 1,b(j) - tmp * b(i));
	}
    }
    for (int i = n; i >= 1; --i) {
	tmp = b(i);
	for (int j = i + 1; j <= n; ++j) tmp -= a(i)(j)*x(j);
	x.set_coord(i - 1,tmp/a(i)(i));
    }
    x.print();
};

template <class T,int n>
void GSMethod(Vector<Vector<T,n>,n> a,Vector<T,n> b,T epsilon) {
    Vector<T,n> x;
    Vector<T,n> y;
    T sum;
    do {
	y = x;
	for (int i = 1; i <= n; ++i) {
	    sum = 0;
	    for (int j = 1; j < i; ++j) sum += -a(i)(j)*x(j);
	    for (int j = i + 1; j <= n; ++j) sum += -a(i)(j)*x(j);
	    sum += b(i);
	    x.set_coord(i-1,sum/a(i)(i));
	}
    } while (fabs(x(1)-y(1)) >= epsilon);

    x.print();
};

#else
// DO NOTHING.
#endif


